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Abstract. We analyze the zero-temperature behavior of the XY Edwards- 
Anderson spin glass model on a square lattice. A newly developed algorithm 
combining exact ground-state computations for Ising variables embedded into the 
planar spins with a specially tailored evolutionary method, resulting in the genetic 
embedded matching (GEM) approach, allows for the computation of numerically 
exact ground states for relatively large systems. This enables a thorough re- 
investigation of the long-standing questions of (i) extensive degeneracy of the 
ground state and (ii) a possible decoupling of spin and chiral degrees of freedom 
in such systems. The new algorithm together with appropriate choices for the 
considered sets of boundary conditions and finite-size scaling techniques allows 
for a consistent determination of the spin and chiral stiffness scaling exponents. 
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1. Introduction 

With their rich behavior at low temperatures, spin glasses take a prominent role in 
the large class of magnetic systems with frustration. The most commonly considered 
Hamiltonian is that of the Edwards- Anderson (EA) model [1] , 



with 0(n) spins Si and random, nearest-neighbor couplings J,j. The wealth of 
behavior of these systems is attributed to the random disorder augmenting the 
frustration effects. Unfortunately, it is precisely this quenched disorder that provides 
an exceptional challenge for the application of the various analytical approximation 
methods well known from the treatment of homogeneous systems. Owing to these 
difficulties, most of the advances in the understanding of spin glass systems beyond 
the celebrated mean-field theory of the Sherrington-Kirkpatrick model [2] have been on 
account of ever more sophisticated numerical simulation and optimization techniques 
[1]. For two-dimensional (2D) systems, where for short-range interactions spin glass 
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order is restricted to zero temperature, the investigation of ground states is of 
prominent interest. In general, finding (exact) ground states of spin glass models 
is a computationally hard problem, where the amount of computer time grows 
exponentially with the size of the system [3]. Here, we explore a new avenue to 
advance methods for the so far much less investigated case of systems with continuous 
spins: we introduce a novel approximate optimization algorithm which, for the 2D 
XY spin glass discussed here, allows to find numerically exact ground states with 
good confidence for systems of up to about 30 x 30 spins [4, 5] . 

Generalizing Peierls' argument for the stability of the ordered phase in 
homogeneous systems to situations with quenched disorder, a droplet scaling theory 
for spin glasses has been formulated [6] . Therein, the role of the droplet surface (free) 
energy is taken on by the width J(L) of the distribution of random couplings for 
a real-space renormalization group decimation at length scale L. In the course of 
renormalization, J(L) scales as J(L) ~ L e % defining the spin stiffness exponent 9 S . 
If the system scales to weak coupling, 9 S < 0, spin-glass order is unstable at finite 
temperature and the system is below its lower critical dimension. This is the situation 
for the EA model in 2D [1], where then 9 S describes the properties of the critical 
point at temperature T = 0, relating the correlation length exponent v = —1/9 S 
[6]. Numerically, the domain- wall free energy might be determined from the energy 
difference between ground states of systems with different types of boundary conditions 
(BCs) chosen such as to induce a relative domain wall [6]. For the n = 1 Ising 
spin glass, the ground-state problem on planar graphs is an exception to the rule, 
being polynomial computationally [3]. Hence, large systems can be treated, leading to 
reliable estimates of 9 S = —0.282(2) (Gaussian J^) resp. 9 S = (bimodal Jy) [1, 7]. 
Due to the presence of strong finite-size corrections, relatively large system sizes and/or 
elaborate finite-size scaling techniques appeared mandatory for consistent estimates 
of 9 S [7, 8]. However, for the case n > 1 of continuous spins, which is more relevant to 
real materials, the lack of effective and efficient algorithms for finding exact ground 
states and the necessary restriction to small systems with L < 12 have led to rather 
less consistent estimates, moving in the range 9 S e [—1, —0.75] [9-11]. 

Moreover, the increased symmetry of the order parameter in the continuous 
spin case has led to speculations about a decoupling of spin and chiral variables 
[12]: since the pattern of frozen spins in the glassy phase has internal degrees of 
freedom, there is a factual difference between proper and improper rotations expressed 
in the decomposition 0(n) = SO(n) x Z2 [13]. The additional Ising like chirality 
variables might order independently of the spins (for systems above their lower critical 
dimension) or, at least, show a different stiffness against fluctuations, resulting in a 
scaling exponent 9 C possibly distinct from 9 S . Indeed, measurements of the chiral 
stiffness for small systems yielded 9 C « —0.38 [10, 11], different from 9 S above. More 
recently, however, Kosterlitz and Akino [14] argued that the choice of BCs in previous 
studies was flawed and they suggest a possibly more appropriate approach leading 
to 9 S w —0.38 « 9 C again for sizes L < 10. The hardly compatible previous results 
for this system hence raise several methodological questions: have numerically exact 
ground states been found?, are the apparent strong finite-size effects under control?, 
have the considered sets of BCs been chosen such as to properly select the intended 
excitations?, and what is the role of scaling corrections explicitly depending on these 
BCs? 
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2. Genetic embedded matching approach 

The treatment of large samples for the 2D Ising case is enabled by a transformation to 
an equivalent problem on the complete graph of frustrated plaquettes: following their 
definition, for each spin configuration frustrated plaquettes have an odd number of 
broken bonds around their perimeter, whereas unfrustrated plaquettes have an even 
number of broken bonds. Thus, drawing "energy strings" dual to the broken bonds, 
these connect pairs of frustrated plaquettes, and the total energy of (1) is (up to 
a constant) given by the total length of energy strings, such that the ground state 
corresponds to a minimum-weight perfect matching of frustrated plaquettes [3]. The 
matching problem can be solved in polynomial time using Edmonds' algorithm [15], 
and for the case of planar graphs its solution is guaranteed to transform back to a 
valid spin configuration [16]. This does not directly apply to the continuous spins 
considered here. We suggest, however, to embed Ising variables into the planar spins 
by decomposing Si = Sf + S^~ = (Si ■ r)r + S-- relative to a random direction r in 
spin space. With respect to reflections of spins along the direction r, the Hamiltonian 
(1) decomposes asW = TT'H + H r ^ with TT'H = - £ (iij> JT. e[ej, and 

jr. = J^Si-rWSj-rl, e\ = sign(S 4 • r). (2) 

Hence, for any fixed r and restricting the movement of spins to reflections along 
r, the Hamiltonian (1) for arbitrary n > 1 takes on the form of an Ising model. 
Consequently, Edmonds' algorithm can be applied to find (one of) the ground state(s) 
of the embedded Ising model. It is obvious that this can never increase the energy 
of the full Hamiltonian (1), but the state found depends on the choice of random 
direction r. To statistically recover the full O(n) symmetry of the Hamiltonian, a 
series of subsequent minimizations is performed with respect to successive random 
choices of r, thus gradually decreasing the total energy via non-local moves. We refer 
to this approach as embedded matching technique. 

It can be shown that, although when the full Hamiltonian (1) is in a ground 
state, all embedded Ising Hamiltonians 7Y r '" must be in one of their respective ground 
states as well, successive minimizations with respect to random directions r are not 
guaranteed to drive the system towards its absolute energy minimum. In other words, 
the non-local embedded matching dynamics has metastable states, but many less than 
the simple local spin quench dynamics used before [9, 11]. To converge to ground states 
with high probability, we insert the embedded matching technique as minimization 
component ("subroutine") in a genetic algorithm [17]: we consider a whole population 
of candidate ground-state configurations and simulate an evolutionary development by 
re-combining (or crossing over) neighboring pairs of parent configurations followed by 
minimization runs for the resulting offspring and replacement of the parents in case of 
lower energy of the offspring. In analogy with the approach of Ref. [17], the crossover 
is performed in a "triadic" fashion, comparing the overlaps with a third, reference 
configuration. This layout is complemented by intermittent random mutation steps 
and performance-guided halving of the population at certain stages to find an optimum 
balance between "genetic" diversity and efficiency of optimization [17]. The choice of 
operation for the crossover of configurations is found to be crucial for the efficiency 
of the approach: it turns out that a random exchange of suitably defined rigid spin 
clusters is far more efficient than an exchange of single spins. These clusters denote 
regions which are highly optimized inside for all configurations of the population (i.e., 
metastable states), but which have to undergo a series of independent rigid O(n) 
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Figure 1. Local rotation matrices between the ground states for a single 20 X 20 
disorder realization with open boundaries relative to domain-wall BCs for spin 
(left) and chiral (right) excitations. 



transformations to make up a true ground state of the system. Careful choice of 
the parameters of the resulting genetic embedded matching (GEM) algorithm and 
application of various statistical tests ensure that indeed independent runs for a single 
given realization of the disorder variables always converge to a state of the same 
energy, up to unprecedented machine precision, which in this way can be guaranteed 
to be a ground state with high reliability [4, 5]. 

We here concentrate on the symmetric, bimodal ± J distribution. For this case 
we find that the ground states computed in independent runs for a single disorder 
realization are always identical to each other up to a global O(n) transformation, 
indicating the lack of accidental degeneracies in contrast to what is found for the 
bimodal Ising case. Hence, after averaging over disorder, the ground state is 
ordered and the ground-state spin correlation function is constant, implying rj = 0. 
To determine the asymptotic ground-state energy per spin e^, ground-states were 
computed for L x L square-lattice systems with L — 6, 8, 10, 12, 16, 20, 24, and 28 for 
open and open-periodic BCs and 5 000 disorder replica per size. Finite-size corrections 
are expected to be purely analytic for the case of open BCs [18] , and a fit to the ansatz 
e(L) = eoo + a/L + b/tf + c/L 3 yields = -1.5520(14) with quality-of-fit Q = 0.35. 
For the open-periodic case, an additional non-analytic term oc is expected to 

occur [18], and a fit of the corresponding data to the form e(L) = e^ + a/L + b/L 2 ~ e 
gives eoo = —1.5525(13), 9 = —0.49(69), Q = 0.35, perfectly consistent with the 
open-boundary result for eoo and, due to the large statistical error, only in qualitative 
agreement with the expected value for the spin-stiffness exponent 6. The resulting eoo 
is about 10% lower than the value eoo = —1.402 of the bimodal Ising spin-glass [18]. 

3. Spin and chiral stiffness exponents 

Conventionally, domain- wall energies have been measured by comparing ground states 
for periodic and antiperiodic (P/AP) BCs [9-11]. In Ref. [14] it was argued, however, 
that the periodicity in both types of BCs forces domain walls into the system, such 
that the corresponding energy difference might not properly capture the energy of 
a single domain wall. There, an improvement is suggested by optimizing over an 




Figure 2. Aspect-ratio scaling of the stiffness exponents 9 S and 9 C for aspect 
ratios R = 1, 2 and 6 as a function of f/il. 



additional global twist variable along the boundary under consideration. Here, to 
start with the cleanest possible setup, in addition to the conventional P/AP BC set, 
we consider open and domain-wall (O/DW) BCs, where for the latter the relative 
orientations of spins linked across the boundary are either tilted by an angle 7r for 
spin domain walls or reflected along an arbitrary but fixed axis for chiral domain walls 
by the introduction of very strong bonds [5, 7]. In Figure 1 we show snapshots of 
spin and chiral excitations forced into the system by the O/DW BCs. To this end, we 
computed a locally averaged 0(2) rotation matrix relating the configurations with O 
and DW BCs and translated it back into a rotation angle (the arrows) and the sign 
of the determinant (—1 for the blue squares). It is apparent that in contrast to the 
discrete Ising case the spin domain walls are rather smeared out and that to a certain 
extent the system appears to relax the spin excitation also through the chiral mode if 
it is found to be softer locally (and vice versa for the chiral excitation). 

From the scaling of the domain-wall energies, [|AS|]j ~ L e , we find a strong 
crossover for the P/AP data from 9 S = -0.724(21) for L < 12 to 6 S = -0.433(26) 
for L > 16, indicating large finite-size effects and a movement from the value found 
for small P/AP computations in previous works [9-11] towards the "optimum twist" 
value of Ref. [14]. The O/DW BCs, on the other hand, yield 6 a = -0.207(12) for spin 
excitations resp. 8 C = —0.090(23) for the chiral domain walls. Hence, although it is 
already clear that the true stiffness exponents are much less negative than estimated 
before, there is still a sizable difference between the P/AP and O/DW results for 9 S , 
indicating incomplete control over finite-size effects. To improve on this, we take into 
account that, due to the independence of BCs for systems in one dimension, corrections 
depending on BCs should disappear as more and more elongated systems are being 
considered [8]. Thus, we additionally performed computations for Lx M systems (the 
change of BCs happening along the edges of length L) with aspect ratios R = M/L = 2 
and 6 with the same statistics. The results are presented in Figure 2 for the case of 
P/AP and O/DW BCs, respectively. Guided by the experience from the Ising case, 
we expect corrections depending on BCs to disappear as 9(R) = 9(R = oo) + Ar/R 
for large R, and indeed the P/AP and O/DW spin data seem to converge for large 
R, a fit to the given form yielding 9 S (R = oo) = -0.338(20) for P/AP BCs and 
9 s (R = oo) = -0.308(30) for O/DW BCs. The O/DW chiral data, on the other hand, 
give C (R = oo) = —0.114(16), clearly distinct from 9 S . 
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4. Conclusions 

Employing a novel "genetic embedded matching algorithm" , we computed numerically 
exact ground states for the 2D XY EA spin glass with ±J couplings and up to 
28 x 28 spins. No accidental degeneracies occur, implying n = 0. Analyzes of the 
domain-wall energies are hampered by strong finite-size effects which, however, can 
be controlled using the aspect-ratio scaling technique. We find consistent estimates of 
9 S = —0.308(30) from different sets of BCs, clearly smaller in modulus than previous 
estimates [9-11, 14], and rather close to 9 S = —0.28 found for the Gaussian 2D Ising 
case. The chiral exponent 6 C = —0.114(16), on the other hand, is found to be clearly 
different from 6 S and closer to 6 S — found for the bimodal 2D Ising spin glass. Note 
also, that our results are rather far from 8 S = —l/v s = —1.0 resp. 6 C = —l/v c = —0.5 
estimated by finite-temperature Monte Carlo simulations [19], which probably is due 
to equilibration problems at low temperatures. 
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